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AN  EFFICIENT  NUMERICAL  MODEL  OF  THE 
PLANETARY  ATMOSPHERIC  BOUNDARY  LAYER 

1.  INTRODUCTION 

There  are  two  general  approaches  to  develop  numerical  models  of 
the  atmospheric  boundary  layer.  The  most  common  approach  has  been  the 
construction  of  a multi-level  boundary  layer  model  (e.g.  Estoque,  1963, 
Sasamori,  1971),  when  detailed  structures  of  the  planetary  boundary 
layer  (PBL)  wind  and  temperature  are  required.  For  global  circulation 
studies,  however,  the  multi-level  PBL  models  may  not  be  practical  in 
view  of  the  computational  demands.  The  alternative  is  to  develop  more 
efficient  models  with  one  or  two  levels  in  the  vertical  specifically 
designed  for  parameterizing  the  bulk  properties  of  the  PBL  for  use  in 
large  scale  atmospheric  circulations,  (Randall,  1976;  Benoit,  1976; 
Stull,  1977).  This  paper  will  describe  an  efficient  two-dimensional, 
time-dependent,  one-level  boundary  layer  model  suitable  for  global 
circulation  studies. 

Two  fundamental  problems  associated  with  the  PBL  parameterization 
are  determinations  of  the  surface  fluxes  and  the  effect  of  flux  pro- 
files in  the  interior  of  the  PBL.  Several  different  schemes  are 
being  used  for  parameterizing  the  boundary  layer  fluxes  (see  e.g. 
Bhumralkar,  1976).  The  simplest  ones  employ  the  usual  bulk  transfer 

relations  with  all  the  transfer  coefficients  for  drag  C and  heat  C 

U n 

assumed  equal  and  prescribed  a prfori.  In  some  cases  different  values 
are  assigned  for  land  and  ocean  surfaces  and  also  some  allowance  is 
made  for  different  stability  conditions.  But  by  and  large,  these 
parameterization  schemes  are  very  crude. 

Note:  Manuscript  submitted  August  30,  1977. 
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Setter  schemes  have  been  formulated  from  similarity  considerations 
of  the  boundary  layer.  The  basic  assumptions  underlying  all  similarity 
theories  are  that  the  boundary  layer  flow  is  horizontally  homogeneous 
and  quasi-stationary,  a very  restrictive  assumption  for  many  real 
atmospheric  situations.  As  pointed  out  by  Arya  (1977),  in  a GCM,  how- 
ever, the  variables  are  considered  to  be  averaged  horizontally  over  a 
fairly  large  grid  area,  and  thus  the  assumption  of  horizontal  homo- 
geneity is  probably  well  justified.  Another  alternative  form  of  the 
similarity  parametric  relation  originally  proposed  by  Deardorff  (1972) 
uses  the  layer-averaged  wind,  temperature  and  humidity.  According  to 

this  approach,  the  drag  and  heat  transfer  coefficients  C-  and  C are 

D H 

calculated  based  on  some  nomograms.  Numerical  results  based  on  this 
approach  will  be  compared  with  those  generated  by  the  Rossby  similarity 
theory. 

Determination  of  the  shapes  of  the  turbulent  flux  profiles  in  the 
interior  of  the  PBL  is  needed  when  the  interface  between  the  GCM  grid 
layers  lie  within  the  PBL.  The  most  common  approach  has  been  to  pre- 
scribe the  flux  profiles  as  a linear  function  of  height.  This  approach 
perhaps  is  adequate  during  the  daytime  convective  hours  when  temperature 
stratifications  are  nearly  adiabatic.  However,  during  other  stability 
conditions,  the  linear  flux  profile  assumption  may  not  be  valid.  Thus, 
another  purpose  of  this  paper  is  to  investigate  the  effect  of  different 
vertical  profiles  of  turbulent  fluxes  on  model  winds  and  temperatures. 

Finally,  we  shall  conclude  this  study  by  testing  the  validity  of 
our  model  by  comparing  it  with  multi-level  boundary  layer  model  results. 
For  this  purpose,  we  shall  simulate  a stably  stratified  air  flow  pass- 
ing from  smooth  to  rough  surfaces. 
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2. 


THE  GOVERNING  EQUATIONS 


The  model  configuration  is  shown  in  Fig.  1.  At  the  top  of  the 
boundary  layer,  i.e.,  z = H,  the  following  boundary  layer  character- 
istics are  assumed:  (1)  the  wind  speed  is  in  geostrophic  balance 
and  the  atmospheric  pressure  is  constant  in  time  and  along  the  x- 
direction.  (2)  The  temperature  and  humidity  fields  remain  unchanged 
throughout  the  whole  integration  period. 

At  the  level  z = H/2,  the  boundary  layer  mean  wind,  temperature, 
and  humidity  are  governed  by  the  following  equations: 
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(-E.)k  « - ^2.  where  k = R/C  . 
dz  p R0  p 


(6) 


All  the  symbols  listed  in  Appendix  I carry  their  usual  physical 
meanings.  In  deriving  this  set  of  equations,  we  have  assumed  that 
the  atmosphere  is  in  hydrostatic  equilibrium. 
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At  the  ground  surface  Z = 0,  the  wind  speed  is  required  to  be 


zero.  The  surface  temperature  T is  determined  by  the  surface  energy 
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budget  equation  which  gives  rise  to  the  following  predictive  equation 
(Bhurmal kar , 1 975) , 


T n+I  T n . 
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Here  = C + (Xc/2u))3  and  C is  the  volumetric  heat  capacity,  X is  the 

thermal  conductivity,  T the  average  daily  surface  temperature,  u>  is 

the  frequency  of  oscillation,  At  the  time  step,  a is  the  Stefan-Bol tzmann 

constant,  C is  the  geostrophic  heat  transfer  coefficient,  C is  the 
H p 

specific  heat  capacity  at  constant  pressure,  and  L is  the  latent  heat 
of  evaporation. 

3.  THE  CLOSURE  PROBLEMS 

A.  PBL  Parameterization  Schemes 

To  close  the  system  of  equations  described  in  Section  2,  the 
flux  terms  in  equations  (1)  - (b)  have  to  be  parameterized.  The  gen- 
eral form  of  the  drag  and  other  transfer  relations  obtained  from 
similarity  arguments  is, 

A 

kU  ,/u  •»  - UnZ  + A) 
h « o 

kV  ^uA  3 - B sign  • f 


k(9.  - 9 )/6.  - - {ZnZ  + C) 
n o o 


k(q  - qJ/q*  ■ * UnZ  + D) 


in  which  A,  B,  C and  0 are  some  similarity  functions  depending  on 
atmospheric  stability.  Recently,  Yamada  (1976)  reexamined  the  similar- 
ity functions  based  on  the  Wangara  data  (Clarke  et  al.,  1971).  Under 
both  the  stable  and  unstable  conditions,  his  proposed  similarity 
functions  show  substantial  improvement  over  the  previous  work  of 
Deardorff  and  Melgarejo  (1975)  and  Arya  (1975).  For  this  study,  we 
shall  use  (8)  to  compute  the  surface  fluxes  based  on  Yamada's  (1976) 
similarity  functions. 

It  should  be  pointed  out  that  the  height  of  the  atmospheric 
boundary  layer  h may  affect  the  similarity  functions.  In  fact,  this 
is  the  factor  differentiating  the  Rossby  number  similarity  theory 
from  the  generalized  similarity  theory.  In  the  Rossby  number  similar- 
aory,  the  boundary  layer  height  is  uniquely  determined  by  a 
ight  Uj./f.  In  the  generalized  similarity  theory,  h is  con- 
,ered  as  an  variable  depending  on  such  factors  as  diurnal  heating, 
large-scale  subsidence,  large-scale  advection  of  heat  and  moisture 
etc..  Recently,  Yu  (1977)  compared  these  two  similarity  theories  and 
concluded  that  for  mid-latitude  boundary  layer  flow,  the  Rossby  number 
similarity  theory  is  just  as  valid  as  in  the  generalized  similarity 
theory.  In  this  paper,  we  adopt  the  Rossby  number  similarity  theory. 

We  shall  compare  our  model  results  with  those  based  on  the  Deardorff 
(1972)  scheme  according  to  which  the  drag  and  heat  transfer  coefficients 

C„  and  C are  calculated  by  some  nomograms  as  shown  in  Fig.  2. 

D H 

B.  Effects  of  PBL  Parameter i zat ion  on  Model  Winds  and 
Temperatures 

In  this  section,  we  shall  test  the  schemes  described  in 
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Section  3A  with  our  simple  model  and  study  their  effects  on  the  model 
winds  and  temperatures.  For  this  purpose,  it  is  sufficient  to  use  the 
one-dimensional  version  of  the  model,  i.e.,  all  the  advective  terms 
in  equations  (1)  - (6)  are  set  to  zero.  The  model  was  integrated  for 
24  hours  with  the  following  initial  and  boundary  conditions: 

Boundary  conditions:  At  Z = H = 500  m 

u=U  = 10  m/s,  v * v =0m/s 
Q 9 

0 = 0H  = 308°  K,  q = qH  = 10  g/kg; 

At  z = 0,  u = v = 0 and  the  surface  temperature  is  predicted  by  (7) ■ 

Initial  conditions:  At  Z = H/2  = 250  m, 
u * 8 m/s,  v * 2m/s 
0 = 305°  K,  q = 12  g/kgj 

At  Z 3 0,  0q  = 300°  K,  qQ  = 14  g/kg 

Numerical  methods  and  finite  difference  equations  are  described 
in  the  Appendix  II.  The  Crank  Nicolson  method  is  used  to  integrate 
the  model  with  a 30  minute  time  step. 

Fig.  3 shows  model  temperatures  calculated  by  the  Rossby  number 
similarity  theory  and  Deardorff  (1972)  schemes  for  two  selected 
heights,  i.e.  Z = 0 and  Z = 2 m.  We  see  that  a difference  as  large 
as  4°  K may  result  when  different  PBL  parameterization  schemes  are 
used.  For  the  model  wind  speed,  a maximum  of  2.5  m/s  difference  is 
found  at  Z = 2 m between  results  calculated  with  the  Rossby  number 
similarity  theory  and  the  Deardorff  (1972)  schemes.  This  is  shown  in 
Fig.  4.  The  difference  in  wind  speed  becomes  smaller,  being  less  than 
I m/s  at  Z = 250  m.  However,  the  effect  on  the  model  temperature  at 
Z = 250  m is  still  significantly  large,  being  on  the  order  of  3-5°  K. 


This  is  illustrated  in  Fig.  5. 

Fig.  6 shows  the  wind  speed  hodograph  at  Z = 250  m for  the  Rossby 
number  similarity  theory  and  Deardorff  (1972)  scheme.  The  veering  of 
the  cross i sobar ic  angles  for  both  cases  are  similar,  but  the  speed  for 
the  Rossby  number  similarity  theory  is  relatively  larger  than  that  for 
the  Deardorff  schemes.  This  may  be  explained  by  the  difference  in 
calculated  surface  fluxes  of  momentum  and  heat  as  shown  in  Table  1. 

We  ran  see  from  Table  1 that  momentum  fluxes  calculated  by  the  Deardorff 
(1972)  scheme  are  much  larger  than  those  calculated  by  the  similarity 
theory.  As  a result,  the  wind  speed  at  Z = 250  m is  much  reduced  with 


3 to  4°  K in  temperature  prediction  occurs  when  different  PBL  para- 
meterization schemes  are  used.  This  points  out  the  need  for  improved 
atmospheric  boundary  layer  parameterization  schemes  in  general 
circulation  models.  Unfortunately  there  is  not  sufficient  data  to 
select  which  scheme  is  best  in  this  study. 

4.  EFFECT  OF  DIFFERENT  VERTICAL  PROFILES  OF  FLUXES  ON 
MODEL  WINDS  AND  TEMPERATURES 

We  shall  assume  the  following  profile  relationships,  i.e., 

C'w'  = C'w’0  exP  {“  m(z/H)2}  (9) 

where  = (u1,  v1 , 0',  q ' ) , the  subscript  o indicates  the  surface 
values,  and  the  parameter  m prescribes  the  shapes  of  exponential 
profiles  as  shown  in  Fig.  7-  From  (9),  it  follows  that  at  Z = H/2, 

£ 1 w 1 = - (m/H)  exp  (-  m/4)  u'w'  . 
d Z 0 

For  1 i near  prof i 1 es , 

- - ^Ts7r  /H 

dZ  O 

Thus,  the  ratio  of  the  vertical  gradient  of  fluxes  between  exponential 
profiles  and  linear  profiles  is  m exp  (-  m/4).  From  Table  2,  it  can 
be  seen  that  changing  the  value  of  m leads  to  the  modification  of 
vertical  flux  gradient  of  the  linear  profile.  For  example,  with  m = 1, 
the  fluxes  calculated  by  (9)  will  be  78  percent  of  that  calculated  by 
linear  profiles,  and  with  m = 4,  the  fluxes  calculated  by  (9)  will  be 
about  147  percent  of  that  calculated  by  linear  profiles.  It  should 
be  noted  that  with  m = 1.4  and  m = 8.6  the  fluxes  calculated  by  (9) 
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l,sMk 


is  nearly  equal  to  that  of  the  linear  profile.  We  integrated  the 
one-dimensional  version  of  the  model  for  24  hours  with  three  different 
values  of  m i.e.,  m = 1,  m * 4 and  m = 8.6  for  the  exponential  flux 
profile  relationship  (9).  It  should  be  noted  that  the  governing 
equations  for  all  the  cases  are  the  same  except  the  flux  divergence 
terms  in  equations  (I)  through  (4)  are  different  due  to  different 
values  of  m. 

Fig.  8 shows  potential  temperatures  at  Z = 250  m calculated  with 
three  different  values  of  m.  We  see  the  difference  in  potential  temp- 
erature between  the  linear  flux  profile  (m  * 8.6)  and  exponential  flux 
profiles  with  m = 1 and  m = 4 is  about  2°  K.  The  difference  in  wind 
speed  is  negligibly  small,  being  less  than  1 m/s  with  the  three 
different  values  of  m.  This  suggests  that  the  shapes  of  turbulent 
flux  profiles  is  not  critical  in  a PBL  parameterization  scheme.  In 
general,  the  use  of  linear  flux  profiles  may  be  sufficient  for  general 
circulation  studies. 

5.  CHANGES  IN  SURFACE  ROUGHNESS 

The  purpose  of  this  section  is  to  determine  if  a simple  model  can 
reproduce  some  of  the  important  features  of  more  sophisticated  multi- 
level boundary  layer  models.  Yu  and  Wagner  (1975)  in  a numerical 
study  concluded  that  when  stably-stratified  air  flow  passes  over  a 
rougher  surface,  the  increased  turbulent  fluxes  of  heat  induced  by  the 
increased  surface  drag  can  Cause  the  formation  of  a heat  island.  The 
model  described  in  Section  2 employs  similar  governing  equations  and 

numerical  schemes  as  those  of  Yu  and  Wagner  (1975).  The  only  differ- 
ence between  these  two  models  is  the  number  of  levels  in  the  vertical. 
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While  our  model  has  only  one  level,  the  model  of  Yu  and  Wagner  (1975) 
had  21  levels.  We  shall  thus  simulate  a stratified  flow  passing  over 
a change  in  surface  roughness  as  shown  in  Table  3-  The  boundary 
conditions  are: 

At  Z * K = 500  m;  u(x,H)  = U =10  m/s; 

9 

v (x , H )=  V = 0;  6(x,  H)  = 0U  = 308°  K; 

g h 

w (x , H)  = 0;  and  p(x,  H )=  P = 900  mb 

rl 

At  z * 0;  u (x , 0)  = v (x , 0)  = w(x,  0)  = 0; 

0(x,  0)  = predicted  by  equation  (7). 

At  the  upwind  lateral  boundary,  the  variables  u,  v,  9 are  assumed 
to  be  uniform  along  the  x-direction,  and  are  obtained  from  the  one- 
dimensional time-dependent  version  of  the  model.  At  the  outflow  lateral 
boundary,  the  upstream  differencing  scheme  does  not  require  boundary 
values.  The  initial  conditions  are: 

At  Z = H/2  = 250  m;  u(x,  z)  = 8 m/s; 
v (x , z)  = 2 m/s;  A(x,  z)  = 305°  K 

o 

At  Z = 0;  9 (x , 0)  = 300  K 

The  finite  difference  equations  are  given  in  Appendix  II.  An 
ir.p’icit  backward-t  ime,  upstream  space  finite  difference  scheme  with 
a 30  minute  time  step  is  used  to  integrate  the  model. 

Fig.  9 shows  the  calculated  diurnal  temperature  variation  at 
two  locations,  i.e.  X = 1 and  X = 5.  Referring  to  Table  3,  we  see  at 
X = 1,  the  roughness  length  is  0.001  m,  while  at  X = 5,  the  roughness 
length  is  2 m.  From  Fig.  9,  it  is  evident  that  over  the  rougher 
surface,  due  to  increased  turbulent  drag  and  heat  flux  transfer,  the 
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air  is  well  mixed  near  the  surface  during  both  day  and  night.  Thus 
the  difference  in  temperature  between  Z = 0 and  Z = 2 m is  very  small 
at  Z * 5.  At  the  upstream  boundary,  i.e.  Z = 1 , temperature  difference 
between  Z * 0 and  Z = 2 m is  large  especially  during  the  nighttime. 

The  most  interesting  feature  is  that  during  the  stable  hours,  (18  - 06) , 
the  temperature  at  X » 5 is  much  warmer  than  that  at  X = 1 . This  is 
attributable  to  the  fact  that  when  a stably-stratified  air  passes  over 
a rougher  surface,  the  energy  is  redistributed  due  to  increased  tur- 
bulence mixing.  This  redistribution  of  heat  results  in  the  formation 
of  a heat  island.  During  the  most  unstable  hours  (12  - lb),  we  notice 
that  the  temperature  near  the  surface  up  to  2 m over  the  rougher  sur- 
face (X  ■ 5)  is  cooler  than  that  at  X = 1 . This  is  due  to  that  inten- 
sive upward  energy  transfer  occuring  over  the  rougher  surface  results 
in  less  energy  available  for  heating  the  air  near  the  surface.  These 
conclusions  in  general  are  in  good  agreement  with  the  numerical  results 
of  Yu  and  Wagner  (1975). 

6.  SUMMARY  AND  CONCLUSIONS 

This  paper  presents  an  efficient  two-dimensional  time  dependent 
one-layer  boundary  layer  model  especially  designed  for  use  in  general 
circulation  models.  The  one-dimensional  version  of  the  model  was 
employed  to  examine  the  effect  of  different  PBL  parameterization 
schemes  on  the  model  winds  and  temperatures.  Two  PBL  parameterization 


schemes  considered  to  be  more  suitable  for  use  in  GCM  studies  were 
tested,  namely,  the  Rossby  number  similarity  theory  and  Deardorff 
(1972)  schemes.  The  effect  of  using  different  PBL  parameterization 
schemes  is  found  to  be  very  significant  on  temperature  prediction  but 


not  so  much  on  the  wind  speed.  Comparison  of  the  numerical  results 
using  these  two  different  PBL  parameterization  schemes  shows  a maximum 
difference  of  A°  K in  the  temperature  prediction.  This  difference  is 
attributable  to  the  large  difference  in  the  determination  of  surface 
momentum  and  heat  fluxes.  The  fluxes  of  momentum  and  heat  calculated 
by  the  similarity  theory  are  much  smaller  than  those  of  Deardorff  (1972) 
scheme.  In  particular,  during  the  unstable  hours,  the  heat  fluxes 
calculated  by  the  Deardorff  (1972)  scheme  are  negative  implying  a upward 
transfer  of  heat  energy.  By  the  similarity  theory,  the  heat  fluxes 
are  calculated  to  be  positive  during  the  unstable  hours  indicating  a 
downward  transfer  of  heat  energy.  This  fundamental  discrepancy  has  to 
be  resolved  by  comparing  results  with  observations.  Clarke  (197*0  did 
a comparison  between  these  two  parameterization  schemes  using  the 
Wangara  <periment  data  and  concluded  that  the  similarity  theory  was 
more  re  le. 

T feet  of  different  vertical  profiles  in  turbulent  fluxes  on 
model  w j and  temperature  was  determined  using  the  one-dimensional 
version  of  the  model.  It  was  concluded  that  the  effect  was  less 
important  and  for  general  circulation  studies,  a linear  profile  might 
be  adequately  sufficient. 

The  two-dimensional  model  was  used  to  simulate  a stably-stratified 
air  flow  passing  over  from  a smooth  to  a rougher  surface.  This  was 
done  to  see  if  our  model  was  capable  of  producing  the  same  physical 
results  as  those  by  more  sophisticated  multi-level  boundary  layer 
models.  It  was  concluded  that  our  numerical  results  agree  essentially 
with  those  of  Yu  and  Wagner  (1975)  in  that  when  stably  stratified 
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flow  passing  over  a rougher  surface,  the  increased  surface  drag  and 
thus  the  increased  turbulent  mixing  can  cause  a heat  island  formation. 
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APPENDIX  I 


List  of  Symbols 


geostrophic  drag  coefficient 

geostrophic  heat  transfer  coefficient 

specific  heat  capacity  at  constant  pressure 

Coriolis  parameter 

gravitational  acceleration 

height  of  the  model  upper  boundary 

height  of  the  model  upper  boundary,  or  sensible  heat  flux 

space  index  for  the  x~direction 

space  index  for  the  z-direction 

time  index  for  the  finite  difference  equations 

atmospheric  pressure 

atmospheric  pressure  of  1000  mb 

specific  humidity  of  air 

specific  humidity  of  air  at  surface 

saturated  specific  humidity  of  air 

surface  friction  humidity 


T surface  temperature 

q 

u,v,w  velocity  components  in  the  x-,  y-,  and  z-d i rect ions , 
respectively 


V vh 


velocity  components  in  the  x-,  y-directions  respectively 
at  Z * h 

surface  friction  velocity 


X,  y,  z 

Z 


coordinates  in  the  downwind,  cross-wind,  and  vertical 
directions  respectively 

surface  roughness  parameter 

normalized  surface  roughness  parameter  («  ZQ/h) 

potential  temperature 

potential  temperature  at  Z = h 

friction  potential  temperature 

air  density 

latest  available  values  in  the  iteration  schemes 
mean  values 

deviations  from  the  mean  values 
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Note  that  using  the  upstream  finite  difference  scheme  requires  no 
outflow  boundary  conditions. 

The  finite  difference  equations  for  u > 0 may  be  written  as, 


Table  I:  Momentum  and  kinematic  heat  fluxes  calculated  by 


Table  2:  The  ratio  between  exponential  flux  profile  and 
linear  flux  profile  for  different  values  of  m 


m 

mexp(-m/4) 

m 

mexp (-m/4) 

0 

0 

5 

1.43 

1 

0.73 

6 

1.33 

1.4 

0.99 

7 

1.22 

1.5 

1.03 

8 

1 .08 

2 

1.21 

8.6 

1.001 

3 

1.42 

9 

0.95 

4 

1.47 

10 

0.88 

Table  3: 

Assumed 

1 ength 

variation  in  x 

-direction  for  rouqhness 

1 

x (km) 

Z (m) 
o 

1 

0 

0.001 

2 

5 

0.001 

3 

10 

1. 00 

4 

15 

1.00 

5 

20 

2.00 

6 

25 

1.00 

7 

30 

1.00 

8 

35 

0.001 

9 

40 

0.001 

10 

45 

0.001 

2 
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- Potential  temperature  at  two  selected  heights,  i.e.,  Z = 0 and  Z»2m  calculated 
by  the  similarity  theory  and  Deardorff  (1972)  scheme 


U-COMPONENT  WIND  SPEED 
AT  z = 250  m 
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MEAN  WIND  SPEED  AT  z = 2m 
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Fig,  4 - Two-meter  mean  wind  speed  and  u-component  wind  speed  at  Z = 250  m 
calculated  by  the  similarity  theory  and  Deardorff  (1972)  scheme 


